# "Media's Influence on LGBTQ Support Across Africa" by Stephen Winkler

# Code to generate:
#   - Figure 2 (First differences of all main models)
#   - Figure A.4a (OLS first differences w/ and w/o social tolerance)
#   - Figure A.4b (Logit first differences w/ and w/o social tolerance)
#   - Figure A.4c (Multilevel first differences w/ and w/o social tolerance)

# For binomial logit & ols first differences, I rerun the models
# with CFE and DCSE. This requires reloading data.
# For the multilelevel models, just need to load the
# results saved as RDS and extract relevant parameters.

rm(list=ls()) # clear environment
source("Rprofile.R") # setup Rprofile (load packages, etc.)

# Load data 
myData <- readRDS("data/afrobarometer.rds") #if using rds file

# Convert Country into Dummy & Rename
dummies <- predict(dummyVars(~ country, data = myData, levelsOnly = TRUE), newdata = myData)
colnames(dummies) <- gsub("country", "", colnames(dummies)) #remove country from column names
colnames(dummies) <- gsub(" ", "", colnames(dummies)) #remove spaces
colnames(dummies) <- gsub("'", "", colnames(dummies)) # remove unusual characters

# Join country dummies back to main dataset
myData <- cbind(myData, dummies)
myData <- dplyr::select(myData, -Algeria, -Sudan, -Egypt) #remove NA countries

# Transform data to numeric b/c sims doesn't like logical/factors
myData$radio <- as.numeric(myData$radio)
myData$tv <- as.numeric(myData$tv)
myData$newspaper <- as.numeric(myData$newspaper)
myData$internet <- as.numeric(myData$internet)
myData$social_media <- as.numeric(myData$social_media)
myData$media_aggregate <- as.numeric(myData$media_aggregate)
myData$media_noradio <- as.numeric(myData$media_noradio)
myData$media_notv <- as.numeric(myData$media_notv)
myData$media_nopaper <- as.numeric(myData$media_nopaper)
myData$media_nointernet <- as.numeric(myData$media_nointernet)
myData$media_nosocialmedia <- as.numeric(myData$media_nosocialmedia)
myData$sexuality2 <- as.numeric(myData$sexuality2)
myData$tolerance_noLGBT2 <- as.numeric(myData$tolerance_noLGBT2)
myData$female <- as.numeric(myData$female)
myData$education <- as.numeric(myData$education)
myData$religiosity <- as.numeric(myData$religiosity)
myData$age <- as.character(myData$age) #age needs to go through character first
myData$age <- as.numeric(myData$age) #then to numeric to maintain real values
myData$income_water <- as.numeric(myData$income_water)
myData$urban <- as.numeric(myData$urban)

#################################
#### Setup all the models   #####
#################################

# Social Media
# Estimate models (OLS & LOGIT)
model_social_media <-
  sexuality2 ~
  social_media +
  media_nosocialmedia +
  tolerance_noLGBT2 + 
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia 

mdata_social_media <- extractdata(model_social_media, myData, na.rm=TRUE, #only extract data used in model
                                  extra = ~respno + district) #add respno & district

logit.social_media <- glm(model_social_media, family = binomial(link="logit"), data = mdata_social_media) #Logit
logit.social_media.robust <- coeftest(logit.social_media, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.social_media <- cluster.vcov(logit.social_media, mdata_social_media$district) #extract clustered vcov matrix
pe.social_media <- logit.social_media.robust[,1] #extract pe from clustered model 

lm.social_media <- lm(model_social_media, data = mdata_social_media) #OLS
lm.social_media.robust <- coeftest(lm.social_media, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.social_media.lm <- cluster.vcov(lm.social_media, mdata_social_media$district) #extract clustered vcov matrix
pe.social_media.lm <- lm.social_media.robust[,1] #extract pe from clustered model 

# Setup model & data (MULTILEVEL)
# Don't have to re-estimate the models b/c results are saved as rds
model.b.social_media <-
  sexuality2 ~
  social_media +
  media_nosocialmedia + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban 
  #(1 + social_media | country) + (1 | district)  #add VSVI for country & VI for District
  # I remove country and district from the model for two reasons
  # 1. it's already estimated so it's not necessary for estimation
  # 2. simcf doesn't like factors, so if I keep them they need to be dummies
mdata_social_media.ml <- extractdata(model.b.social_media, myData, na.rm = TRUE,
                                     extra = ~country) #only extract data used in model

ml.social_media <- readRDS("analysis/results/social_media.b.result.rds") # Multilevel results
pe.social_media.ml <- fixef(ml.social_media) # point estimates
vc.social_media.ml <- vcov(ml.social_media) # var-cov matrix

# Internet
# Estimate models (OLS & LOGIT)
model_internet <-
  sexuality2 ~
  internet +
  media_nointernet +
  tolerance_noLGBT2 +
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata_int <- extractdata(model_internet, myData, na.rm=TRUE, #only extract data used in model
                                  extra = ~respno + district) #add respno & district

logit.internet <- glm(model_internet, family = binomial(link="logit"), data = mdata_int) #Logit
logit.int.robust <- coeftest(logit.internet, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.int <- cluster.vcov(logit.internet, mdata_int$district) #extract clustered vcov matrix
pe.int <- logit.int.robust[,1] #extract pe from clustered model 

lm.internet <- lm(model_internet, data = mdata_int) #OLS
lm.int.robust <- coeftest(lm.internet, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.int.lm <- cluster.vcov(lm.internet, mdata_int$district) #extract clustered vcov matrix
pe.int.lm <- lm.int.robust[,1] #extract pe from clustered model 

# Setup model & data (MULTILEVEL)
# Don't have to re-estimate the models b/c results are saved as rds
model.b.internet <-
  sexuality2 ~
  internet +
  media_nointernet + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban 
  #(1 + internet | country) + (1 | district)  #add VSVI for country & VI for District
mdata_int.ml <- extractdata(model.b.internet, myData, na.rm = TRUE) #only extract data used in model

ml.internet <- readRDS("analysis/results/internet.b.result.rds") # Multilevel results
pe.int.ml <- fixef(ml.internet) # point estimates
vc.int.ml <- vcov(ml.internet) # var-cov matrix

# Newspaper
# Estimate models (OLS & LOGIT)
model_newspaper <-
  sexuality2 ~
  newspaper +
  media_nopaper +
  tolerance_noLGBT2 +
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata_paper <- extractdata(model_newspaper, myData, na.rm=TRUE, #only extract data used in model
                         extra = ~respno + district) #add respno & district

logit.paper <- glm(model_newspaper, family = binomial(link="logit"), data = mdata_paper) #Logit
logit.paper.robust <- coeftest(logit.paper, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.paper <- cluster.vcov(logit.paper, mdata_paper$district) #extract clustered vcov matrix
pe.paper <- logit.paper.robust[,1] #extract pe from clustered model 

lm.paper <- lm(model_newspaper, data = mdata_paper) #OLS
lm.paper.robust <- coeftest(lm.paper, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.paper.lm <- cluster.vcov(lm.paper, mdata_paper$district) #extract clustered vcov matrix
pe.paper.lm <- lm.paper.robust[,1] #extract pe from clustered model 

# Setup model & data (MULTILEVEL)
# Don't have to re-estimate the models b/c results are saved as rds
model.b.newspaper <-
  sexuality2 ~
  newspaper +
  media_nopaper + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban 
  #(1 + newspaper | country) + (1 | district)  #add VSVI for country & VI for District
mdata_paper.ml <- extractdata(model.b.newspaper, myData, na.rm = TRUE) #only extract data used in model

ml.newspaper <- readRDS("analysis/results/newspaper.b.result.rds") # Multilevel results
pe.paper.ml <- fixef(ml.newspaper) # point estimates
vc.paper.ml <- vcov(ml.newspaper) # var-cov matrix

# TV
# Estimate models (OLS & LOGIT)
model_tv <-
  sexuality2 ~
  tv +
  media_notv +
  tolerance_noLGBT2 +
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata_tv <- extractdata(model_tv, myData, na.rm=TRUE, #only extract data used in model
                           extra = ~respno + district) #add respno & district

logit.tv <- glm(model_tv, family = binomial(link="logit"), data = mdata_tv) #Logit
logit.tv.robust <- coeftest(logit.tv, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.tv <- cluster.vcov(logit.tv, mdata_tv$district) #extract clustered vcov matrix
pe.tv <- logit.tv.robust[,1] #extract pe from clustered model 

lm.tv <- lm(model_tv, data = mdata_tv) #OLS
lm.tv.robust <- coeftest(lm.tv, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.tv.lm <- cluster.vcov(lm.tv, mdata_tv$district) #extract clustered vcov matrix
pe.tv.lm <- lm.tv.robust[,1] #extract pe from clustered model 

# Setup model & data (MULTILEVEL)
# Don't have to re-estimate the models b/c results are saved as rds
model.b.tv <-
  sexuality2 ~
  tv +
  media_notv + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban 
  #(1 + tv | country) + (1 | district)  #add VSVI for country & VI for District
mdata_tv.ml <- extractdata(model.b.tv, myData, na.rm = TRUE) #only extract data used in model
ml.tv <- readRDS("analysis/results/tv.b.result.rds") # Multilevel results
pe.tv.ml <- fixef(ml.tv) # point estimates
vc.tv.ml <- vcov(ml.tv) # var-cov matrix

# Radio
# Estimate models (OLS & LOGIT)
model_radio <-
  sexuality2 ~
  radio +
  media_noradio +
  tolerance_noLGBT2 +
  female + education + religiosity + age + income_water + urban +
  Benin + Botswana + BurkinaFaso + Burundi + Cameroon + 
  CapeVerde + CotedIvoire + Gabon + Ghana + Guinea + Kenya + 
  Lesotho + Liberia + Madagascar + Malawi + Mali + Mauritius + Morocco + 
  Mozambique + Namibia + Niger + Nigeria + SãoToméandPríncipe + Senegal +
  SierraLeone + SouthAfrica + Swaziland + Tanzania + Togo + Tunisia +
  Uganda + Zambia

mdata_radio <- extractdata(model_radio, myData, na.rm=TRUE, #only extract data used in model
                        extra = ~respno + district) #add respno & district

logit.radio <- glm(model_radio, family = binomial(link="logit"), data = mdata_radio) #Logit
logit.radio.robust <- coeftest(logit.radio, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.radio <- cluster.vcov(logit.radio, mdata_radio$district) #extract clustered vcov matrix
pe.radio <- logit.radio.robust[,1] #extract pe from clustered model 

lm.radio <- lm(model_radio, data = mdata_radio) #OLS
lm.radio.robust <- coeftest(lm.radio, vcov. = function(x) cluster.vcov(x, ~ district)) #add clust std err
vc.radio.lm <- cluster.vcov(lm.radio, mdata_radio$district) #extract clustered vcov matrix
pe.radio.lm <- lm.radio.robust[,1] #extract pe from clustered model 

# Setup model & data (MULTILEVEL)
# Don't have to re-estimate the models b/c results are saved as rds
model.b.radio <-
  sexuality2 ~
  radio +
  media_noradio + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban 
  #(1 + radio | country) + (1 | district)  #add VSVI for country & VI for District
mdata_radio.ml <- extractdata(model.b.radio, myData, na.rm = TRUE) #only extract data used in model
ml.radio <- readRDS("analysis/results/radio.b.result.rds") # Multilevel results
pe.radio.ml <- fixef(ml.radio) # point estimates
vc.radio.ml <- vcov(ml.radio) # var-cov matrix

#################################
## Create Meaningful
## First Difference Scenarios
## And run simulations 
#################################

# Initialize scenarios to mean values of covariates (OLS & LOGIT)
s_media.scen <- cfMake(model_social_media, data=mdata_social_media, nscen=1)
int.scen <- cfMake(model_internet, data=mdata_int, nscen=1)
paper.scen <- cfMake(model_newspaper, data=mdata_paper, nscen=1)
tv.scen <- cfMake(model_tv, data=mdata_tv, nscen=1)
radio.scen <- cfMake(model_radio, data=mdata_radio, nscen=1)

# Initialize scenarios to mean values of covariates (MULTILEVEL)
s_media.scen.ml <- cfMake(model.b.social_media, data=mdata_social_media.ml, nscen=1)
int.scen.ml <- cfMake(model.b.internet, data=mdata_int.ml, nscen=1)
paper.scen.ml <- cfMake(model.b.newspaper, data=mdata_paper.ml, nscen=1)
tv.scen.ml <- cfMake(model.b.tv, data=mdata_tv.ml, nscen=1)
radio.scen.ml <- cfMake(model.b.radio, data=mdata_radio.ml, nscen=1)

# Configure scenario:  Social Media from 1 to 5
s_media.scen <- cfName(s_media.scen, "social_media 1 vs 5", scen=1)
s_media.scen <- cfChange(s_media.scen, "social_media", 
                     x = 5,
                     xpre=1,
                     scen=1)
# Configure scenario:  Social Media from 1 to 5 (ML)
s_media.scen.ml <- cfName(s_media.scen.ml, "social_media 1 vs 5", scen=1)
s_media.scen.ml <- cfChange(s_media.scen.ml, "social_media", 
                         x = 5,
                         xpre=1,
                         scen=1)
# Configure scenario:  Internet from 1 to 5
int.scen <- cfName(int.scen, "internet 1 vs 5", scen=1)
int.scen <- cfChange(int.scen, "internet", 
                     x = 5,
                     xpre=1,
                     scen=1)
# Configure scenario:  Internet from 1 to 5 (ML)
int.scen.ml <- cfName(int.scen.ml, "internet 1 vs 5", scen=1)
int.scen.ml <- cfChange(int.scen.ml, "internet", 
                     x = 5,
                     xpre=1,
                     scen=1)
# Configure scenario:  Newspaper from 1 to 5
paper.scen <- cfName(paper.scen, "newspaper 1 vs 5", scen=1)
paper.scen <- cfChange(paper.scen, "newspaper", 
                     x = 5,
                     xpre=1,
                     scen=1)
# Configure scenario:  Newspaper from 1 to 5 (ML)
paper.scen.ml <- cfName(paper.scen.ml, "newspaper 1 vs 5", scen=1)
paper.scen.ml <- cfChange(paper.scen.ml, "newspaper", 
                       x = 5,
                       xpre=1,
                       scen=1)
# Configure scenario:  TV from 1 to 5
tv.scen <- cfName(tv.scen, "tv 1 vs 5", scen=1)
tv.scen <- cfChange(tv.scen, "tv", 
                     x = 5,
                     xpre=1,
                     scen=1)
# Configure scenario:  TV from 1 to 5 (ML)
tv.scen.ml <- cfName(tv.scen.ml, "tv 1 vs 5", scen=1)
tv.scen.ml <- cfChange(tv.scen.ml, "tv", 
                    x = 5,
                    xpre=1,
                    scen=1)
# Configure scenario:  Radio from 1 to 5
radio.scen <- cfName(radio.scen, "radio 1 vs 5", scen=1)
radio.scen <- cfChange(radio.scen, "radio", 
                       x = 5,
                       xpre=1,
                       scen=1)
# Configure scenario:  Radio from 1 to 5 (ML)
radio.scen.ml <- cfName(radio.scen.ml, "radio 1 vs 5", scen=1)
radio.scen.ml <- cfChange(radio.scen.ml, "radio", 
                       x = 5,
                       xpre=1,
                       scen=1)

# Simulate conditional expectations for these counterfactuals
sims <- 10000

# Logit, Linear, ML simulations
# Social media
simbetas.s_media <- mvrnorm(sims, pe.social_media, vc.social_media) #logit 
s_media.qoi <- logitsimfd(s_media.scen, simbetas.s_media, ci=0.95)
simbetas.s_media.lm <- mvrnorm(sims, pe.social_media.lm, vc.social_media.lm) #linear
s_media.qoi.lm <- linearsimfd(s_media.scen, simbetas.s_media.lm, ci=0.95)
simbetas.s_media.ml <- mvrnorm(sims, pe.social_media.ml, vc.social_media.ml) #multilevel
s_media.qoi.ml <- logitsimfd(s_media.scen.ml, simbetas.s_media.ml, ci=0.95)

# Internet
simbetas.int <- mvrnorm(sims, pe.int, vc.int) #logit
int.qoi <- logitsimfd(int.scen, simbetas.int, ci=0.95)
simbetas.int.lm <- mvrnorm(sims, pe.int.lm, vc.int.lm) #linear
int.qoi.lm <- linearsimfd(int.scen, simbetas.int.lm, ci=0.95)
simbetas.int.ml <- mvrnorm(sims, pe.int.ml, vc.int.ml) #multilevel
int.qoi.ml <- logitsimfd(int.scen.ml, simbetas.int.ml, ci=0.95)

# Newspaper
simbetas.paper <- mvrnorm(sims, pe.paper, vc.paper) #logit
paper.qoi <- logitsimfd(paper.scen, simbetas.paper, ci=0.95)
simbetas.paper.lm <- mvrnorm(sims, pe.paper.lm, vc.paper.lm) #linear
paper.qoi.lm <- linearsimfd(paper.scen, simbetas.paper.lm, ci=0.95)
simbetas.paper.ml <- mvrnorm(sims, pe.paper.ml, vc.paper.ml) #multilevel
paper.qoi.ml <- logitsimfd(paper.scen.ml, simbetas.paper.ml, ci=0.95)

# TV
simbetas.tv <- mvrnorm(sims, pe.tv, vc.tv) #logit
tv.qoi <- logitsimfd(tv.scen, simbetas.tv, ci=0.95)
simbetas.tv.lm <- mvrnorm(sims, pe.tv.lm, vc.tv.lm) #linear
tv.qoi.lm <- linearsimfd(tv.scen, simbetas.tv.lm, ci=0.95)
simbetas.tv.ml <- mvrnorm(sims, pe.tv.ml, vc.tv.ml) #multilevel
tv.qoi.ml <- logitsimfd(tv.scen.ml, simbetas.tv.ml, ci=0.95)

# Radio
simbetas.radio <- mvrnorm(sims, pe.radio, vc.radio) #logit
radio.qoi <- logitsimfd(radio.scen, simbetas.radio, ci=0.95)
simbetas.radio.lm <- mvrnorm(sims, pe.radio.lm, vc.radio.lm) #linear
radio.qoi.lm <- logitsimfd(radio.scen, simbetas.radio.lm, ci=0.95)
simbetas.radio.ml <- mvrnorm(sims, pe.radio.ml, vc.radio.ml) #multilevel
radio.qoi.ml <- logitsimfd(radio.scen.ml, simbetas.radio.ml, ci=0.95)

# Convert to Dataframe so I can plot w/ ggplot
# IMPORTANT NOTE: Figure A.4a-c in the Appendix plots the results from these models
# with and without a control for social tolerance. To create this figure, you must
# comment out the tolerance variable in every model above and re-run the model. 

# logit fd
s_media.fd.logit <- as.data.frame(s_media.qoi) %>% mutate(medium = "social_media", model = "logit") 
int.fd.logit <- as.data.frame(int.qoi) %>% mutate(medium = "internet", model = "logit")
paper.fd.logit <- as.data.frame(paper.qoi) %>% mutate(medium = "newspaper", model = "logit")
tv.fd.logit <- as.data.frame(tv.qoi) %>% mutate(medium = "tv", model = "logit")
radio.fd.logit <- as.data.frame(radio.qoi) %>% mutate(medium = "radio", model = "logit")

# logit fd w/o tolerance
s_media.fd.logit.notol <- as.data.frame(s_media.qoi) %>% mutate(medium = "social_media", model = "logit_no_tolerance") 
int.fd.logit.notol <- as.data.frame(int.qoi) %>% mutate(medium = "internet", model = "logit_no_tolerance")
paper.fd.logit.notol <- as.data.frame(paper.qoi) %>% mutate(medium = "newspaper", model = "logit_no_tolerance")
tv.fd.logit.notol <- as.data.frame(tv.qoi) %>% mutate(medium = "tv", model = "logit_no_tolerance")
radio.fd.logit.notol <- as.data.frame(radio.qoi) %>% mutate(medium = "radio", model = "logit_no_tolerance")

# linear fd
s_media.fd.lm <- as.data.frame(s_media.qoi.lm) %>% mutate(medium = "social_media", model = "ols")
int.fd.lm <- as.data.frame(int.qoi.lm) %>% mutate(medium = "internet", model = "ols")
paper.fd.lm <- as.data.frame(paper.qoi.lm) %>% mutate(medium = "newspaper", model = "ols")
tv.fd.lm <- as.data.frame(tv.qoi.lm) %>% mutate(medium = "tv", model = "ols")
radio.fd.lm <- as.data.frame(radio.qoi.lm) %>% mutate(medium = "radio", model = "ols")

# linear fd w/o tolerance
s_media.fd.lm.notol <- as.data.frame(s_media.qoi.lm) %>% mutate(medium = "social_media", model = "ols_no_tolerance")
int.fd.lm.notol <- as.data.frame(int.qoi.lm) %>% mutate(medium = "internet", model = "ols_no_tolerance")
paper.fd.lm.notol <- as.data.frame(paper.qoi.lm) %>% mutate(medium = "newspaper", model = "ols_no_tolerance")
tv.fd.lm.notol <- as.data.frame(tv.qoi.lm) %>% mutate(medium = "tv", model = "ols_no_tolerance")
radio.fd.lm.notol <- as.data.frame(radio.qoi.lm) %>% mutate(medium = "radio", model = "ols_no_tolerance")

# multilevel fd
s_media.fd.ml <- as.data.frame(s_media.qoi.ml) %>% mutate(medium = "social_media", model = "multilevel")
int.fd.ml <- as.data.frame(int.qoi.ml) %>% mutate(medium = "internet", model = "multilevel")
paper.fd.ml <- as.data.frame(paper.qoi.ml) %>% mutate(medium = "newspaper", model = "multilevel")
tv.fd.ml <- as.data.frame(tv.qoi.ml) %>% mutate(medium = "tv", model = "multilevel")
radio.fd.ml <- as.data.frame(radio.qoi.ml) %>% mutate(medium = "radio", model = "multilevel")

# multilevel fd w/o tolerance
s_media.fd.ml.notol <- as.data.frame(s_media.qoi.ml) %>% mutate(medium = "social_media", model = "multilevel_no_tolerance")
int.fd.ml.notol <- as.data.frame(int.qoi.ml) %>% mutate(medium = "internet", model = "multilevel_no_tolerance")
paper.fd.ml.notol <- as.data.frame(paper.qoi.ml) %>% mutate(medium = "newspaper", model = "multilevel_no_tolerance")
tv.fd.ml.notol <- as.data.frame(tv.qoi.ml) %>% mutate(medium = "tv", model = "multilevel_no_tolerance")
radio.fd.ml.notol <- as.data.frame(radio.qoi.ml) %>% mutate(medium = "radio", model = "multilevel_no_tolerance")

fd_logit <- bind_rows(s_media.fd.logit, int.fd.logit, paper.fd.logit,
                                   tv.fd.logit, radio.fd.logit)

fd_lm <- bind_rows(s_media.fd.lm, int.fd.lm, paper.fd.lm,
                   tv.fd.lm, radio.fd.lm)

fd_ml <- bind_rows(s_media.fd.ml, int.fd.ml, paper.fd.ml,
                      tv.fd.ml, radio.fd.ml)

fd_all <- bind_rows(fd_logit, fd_lm, fd_ml)

####
fd_logit_no_tolerance <- bind_rows(s_media.fd.logit.notol, int.fd.logit.notol, paper.fd.logit.notol,
                                   tv.fd.logit.notol, radio.fd.logit.notol)

fd_lm_no_tolerance <- bind_rows(s_media.fd.lm.notol, int.fd.lm.notol, paper.fd.lm.notol,
                   tv.fd.lm.notol, radio.fd.lm.notol)

fd_ml_no_tolerance <- bind_rows(s_media.fd.ml.notol, int.fd.ml.notol, paper.fd.ml.notol,
                                tv.fd.ml.notol, radio.fd.ml.notol)

all_logit <- bind_rows(fd_logit, fd_logit_no_tolerance)

all_lm <- bind_rows(fd_lm, fd_lm_no_tolerance)

all_ml <- bind_rows(fd_ml, fd_ml_no_tolerance)

#################
## CREATE PLOTS
#################

# Just Logit w/ tolerance 
# set levels
fd_logit$medium <- factor(fd_logit$medium,
                        levels = c("radio",
                                   "tv",
                                   "newspaper",
                                   "internet",
                                   "social_media")) 

fd_plot_logit <- ggplot(fd_logit) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = medium
                      #color = (model),
                      #group = (model),
                      #shape = (model)
                      )) + 
  labs(x = "", y = "difference in probability of supporting LGBTQs") +
  ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' consumption") + 
  coord_flip() + 
  scale_x_discrete(labels=c("Radio", "TV", "Newspaper", "Internet", "Social media")) + 
  scale_y_continuous(labels = scales::percent, breaks = c(-0.01, 0.0, 0.02, 0.04)) +
  #scale_color_manual(name = "Model",
  #                   values = (brewer.pal(4, "Accent"))) +
  #scale_shape_manual(name = "Model",
  #                   values = rev(c(0,1,2))) +
  #ylim(-0.015,0.06, labels=scales::percent) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5),
                     legend.position = c(.9,.5))
fd_plot_logit

ggsave(filename = "figures/fd_plot_logit.pdf", #Save it to directory
       plot = fd_plot,
       width = 8, height = 6)

# Create plot with Logit, ols, multilevel w/ tolerance
# Figure 2 in Main Text 

# set levels
fd_all$medium <- factor(fd_all$medium,
                        levels = c("radio",
                                   "tv",
                                   "newspaper",
                                   "internet",
                                   "social_media")) 
fd_all %<>% mutate(model = as.factor(model),
                   model = fct_relevel(model,
                               c("multilevel", "logit", "ols")) )

dodge <- position_dodge(width =.3) #set offset positioning
fd_plot_all <- ggplot(fd_all) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = medium,
                      color = (model),
                      #group = (model),
                      shape = (model)),
                  position = dodge) + 
  labs(x = "", y = "difference in probability of supporting LGBTQs (95% CI)") +
  #ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' consumption") + 
  coord_flip() + 
  scale_x_discrete(labels=c("Radio", "TV", "Newspaper", "Internet", "Social media")) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 2), breaks = c(-0.02, 0.0, 0.02, 0.04, 0.06)) +
  scale_color_manual(name = "Model",
                     values = c("gray80", "gray50", "black")) +
  scale_shape_manual(name = "Model",
                     values = rev(c(0,1,2))) +
  #ylim(-0.015,0.06, labels=scales::percent) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5),
                     legend.position = c(.9,.5))
fd_plot_all

ggsave(filename = "figures/fd_plot_all.pdf", #Save it to directory
       plot = fd_plot_all,
       width = 8, height = 6)


# Logit w/ and w/o tolerance
# Figure A.4b in Appendix 

# set levels
all_logit$medium <- factor(all_logit$medium,
                        levels = c("radio",
                                   "tv",
                                   "newspaper",
                                   "internet",
                                   "social_media")) 

dodge <- position_dodge(width =.3) #set offset positioning
fd_plot_all_logit <- ggplot(all_logit) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = medium,
                      color = (model),
                      #group = (model),
                      shape = (model)),
                  position = dodge) + 
  labs(x = "", y = "difference in probability of supporting LGBTQs") +
  ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' consumption") + 
  coord_flip() + 
  scale_x_discrete(labels=c("Radio", "TV", "Newspaper", "Internet", "Social media")) + 
  scale_y_continuous(labels = scales::percent, breaks = c(-0.01, 0.0, 0.02, 0.04)) +
  scale_color_manual(name = "Model",
                     values = c("gray50", "black")) +
  scale_shape_manual(name = "Model",
                     values = rev(c(0,1))) +
  #ylim(-0.015,0.06, labels=scales::percent) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5),
                     legend.position = c(.9,.5))
fd_plot_all_logit

ggsave(filename = "figures/fd_plot_logit_comparison.pdf", #Save it to directory
       plot = fd_plot_all_logit,
       width = 8, height = 6)

# OLS w/ and w/o tolerance
# Figure A.4a in Appendix 

# set levels
all_lm$medium <- factor(all_lm$medium,
                           levels = c("radio",
                                      "tv",
                                      "newspaper",
                                      "internet",
                                      "social_media")) 

dodge <- position_dodge(width =.3) #set offset positioning
fd_plot_all_lm <- ggplot(all_lm) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = medium,
                      color = (model),
                      #group = (model),
                      shape = (model)),
                  position = dodge) + 
  labs(x = "", y = "difference in probability of supporting LGBTQs") +
  ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' consumption") + 
  coord_flip() + 
  scale_x_discrete(labels=c("Radio", "TV", "Newspaper", "Internet", "Social media")) + 
  scale_y_continuous(labels = scales::percent, breaks = c(-0.01, 0.0, 0.02, 0.04)) +
  scale_color_manual(name = "Model",
                     values = c("gray50", "black")) +
  scale_shape_manual(name = "Model",
                     values = rev(c(0,1))) +
  #ylim(-0.015,0.06, labels=scales::percent) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5),
                     legend.position = c(.9,.5))
fd_plot_all_lm

ggsave(filename = "figures/fd_plot_lm_comparison.pdf", #Save it to directory
       plot = fd_plot_all_lm,
       width = 8, height = 6)

# Multilevel w/ and w/o tolerance
# Figure A.4c in Appendix

# set levels
all_ml$medium <- factor(all_ml$medium,
                        levels = c("radio",
                                   "tv",
                                   "newspaper",
                                   "internet",
                                   "social_media")) 

dodge <- position_dodge(width =.3) #set offset positioning
fd_plot_all_ml <- ggplot(all_ml) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = medium,
                      color = (model),
                      #group = (model),
                      shape = (model)),
                  position = dodge) + 
  labs(x = "", y = "difference in probability of supporting LGBTQs") +
  ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' consumption") + 
  coord_flip() + 
  scale_x_discrete(labels=c("Radio", "TV", "Newspaper", "Internet", "Social media")) + 
  scale_y_continuous(labels = scales::percent, breaks = c(-0.01, 0.0, 0.02, 0.04)) +
  scale_color_manual(name = "Model",
                     values = c("gray50", "black")) +
  scale_shape_manual(name = "Model",
                     values = rev(c(0,1))) +
  #ylim(-0.015,0.06, labels=scales::percent) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5),
                     legend.position = c(.9,.5))
fd_plot_all_ml

ggsave(filename = "figures/fd_plot_ml_comparison.pdf", #Save it to directory
       plot = fd_plot_all_ml,
       width = 8, height = 6)
